Accurate black hole evolutions by fourth-order numerical relativity 
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We present techniques for successfully performing numerical relativity simulations of binary black 
holes with fourth-order accuracy. Our simulations are based on a new coding framework which cur- 
rently supports higher order finite differencing for the BSSN formulation of Einstein's equations, but 
which is designed to be readily applicable to a broad class of formulations. We apply our techniques 
to a standard set of numerical relativity test problems, demonstrating the fourth-order accuracy of 
the solutions. Finally we apply our approach to binary black hole head-on collisions, calculating 
the waveforms of gravitational radiation generated and demonstrating significant improvements in 
waveform accuracy over second-order methods with typically achievable numerical resolution. 
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I. INTRODUCTION 

Gravitational wave astronomy will soon provide astro- 
physicists with a new tool to observe and analyze some 
of the most energetic phenomena in our universe. Colli- 
sions of compact objects, such as neutron stars, stellar- 
mass black holes, and super-massive black holes, should 
produce characteristic gravitational wave signals that are 
observable to cosmological distances. In particular bi- 
nary black hole systems are among the most promising 
sources of gravitational waves for both the current gen- 
eration of ground-based detectors, such as LIGO and 
for the next generation of space-based detectors, such as 
LISA While early ground-based detectors will need 

to use templates from theoretically derived waveforms in 
order to extract the weak signal from the much larger 
noise, space-based detectors will require accurate models 
of the gravitational wave sources to extract the charac- 
teristic information (e.g., mass, spin) from the detected 
gravitational waves. Thus there is a critical need for ac- 
curate gravitational waveforms from realistic simulations 
of merging black holes. 

Numerical studies of Einstein's equations play an es- 
sential role in studying highly dynamic systems, such 
as binary black holes, that have little or no symme- 
try. These are not only computationally very demanding, 
requiring high-performance massive parallel computers, 
but also mathematically and numerically very challeng- 
ing. Despite these difficulties a great deal of progress 
has been made in the last few years p|, and it is now 
possible to follow binary black hole evolutions during the 
last moments of the merger phase H 1 and perhaps 
even through a complete orbit The calculation of 
the gravitational radiation emission from merging binary 
black holes has also been made possible through the use 
of the Lazarus approach, which bridges numerical relativ- 
ity and perturbative techniques to extract approximate 
gravitational waveforms [9J, llOl llll \VA Il3j . These cal- 
culations are in good agreement with recent numerical 
relativity evolutions 0. 



At present, one of the most serious limiting factor in 
numerical relativity has been the inability to extract ac- 
curate and reliable results from stable three-dimensional 
numerical simulations of coalescing binary black hole 
spacetimes. Considerable resources have been dedicated 
to finding numerically stable formulations of the Einstein 
evolution equations HQHHElEHIiillil 
|H |H Hit This effort has led t0 stable evolutions of 
single black hole spacetimes [2^, |2(| |27l H^. However, 
there is a growing realization that solving a well-posed 
formulation of the Einstein equations with second-order 
accurate finite differencing is not sufficient to produce ac- 
curate simulations of binary black hole spacetimes. Un- 
fortunately, black hole simulations require large domains 
with high resolution near the horizons, and the unigrid 
second-order accurate codes developed in the past are 
not sufficient, given the foreseeable constraints on mem- 
ory and CPU speed, to p roduce accurate waveforms from 
merging black holes |29| |. The Numerical Relativity com- 
munity is currently pursuing several different approaches 
to produce accurate evolutions. Among them are the 
use of finite elements, spectral methods JEIMDj' adap- 
tive mesh finite-difference codes [32, l33t l34l l35l l36j | , and 
higher order finite-difference codes j33|. In addition, 
fourth-order accurate algorithms have recently been de- 
signed to evolve perturbations of nonrotating |38j and 
rotating black holes 39] . Ideally, the next generation of 
finite-difference codes should combine higher order finite- 
differencing in combination with mesh refinement tech- 
niques. However, considerable progress toward highly ac- 
curate evolutions can be achieved using a unigrid higher- 
order finite-difference code along with a coordinate sys- 
tem, such as 'Fisheye' [T^j, that concentrates the grid- 
points in the central region containing the black holes. In 
this paper we present a successful approach to performing 
fourth-order accurate numerical relativity simulations. 

Our simulations are based on a new numerical rela- 
tivity evolution code framework, LazEv, to solve the full 
non-linear Einstein's equations in 3D. The LazEv frame- 
work is designed to be modular. The code consists of 
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a 'Method of Lines' time integrator along with several 
Mathematica scripts |4(j that convert tensorial equations 
into finite-difference algorithms (of arbitrary order) in C. 
This modularity allows us to quickly implement new for- 
mulations and gauge choices as needed. We have cur- 
rently implemented the ADM formulation as well as sev- 
eral 'flavors' of the BSSN formulation [3 E E3- In 
Sec. [n] we describe the LazEv framework in detail. In 
Sec. IIIII we enumerate many of the basic techniques uti- 
lized in numerical relativity simulations which we have 
implemented within LazEv, providing the foundation for 
detailed testing and application in black hole simulations. 

In Scc. lIVI wc test our fourth-order techniques on stan- 
dard numerical relativity testbed problems demon- 
strating fourth-order accuracy. Finally, in Sec we ap- 
ply these methods to study head-on binary black hole 
collisions. Although head-on collisions have limited as- 
trophysical relevance, there has been recent interest in us- 
ing these spacetimes to develop and test numerical tech- 
niques for evolving more generic binary black hole space- 
times [33l 03 • In this paper we demonstrate that our 
fourth-order techniques can realize significant improve- 
ments in gravitational radiation waveform accuracies at 
typically achievable numerical resolutions. 



II. THE LAZEV EVOLUTION FRAMEWORK 

The LazEv evolution framework is a foundation for 
rapid implementation of numerical relativity evolution 
system routines, supporting higher order finite differ- 
encing techniques with 'method of lines' (MoL) integra- 
tion j43|. Our numerical code uses the Cactus Compu- 
tational Toolkit [44| for parallelization and 10. We have 
constructed the code so that it is compatible with the ini- 
tial data and analysis thorns provided with Cactus. Spe- 
cific evolution routines are produced using Mathematica 
scripts ^3 that convert tensorial equations into finite 
difference algorithms in C. These scripts provide several 
types of consistency checks that prevent many types of 
potential coding errors in setting up the tensorial equa- 
tions. The Mathematica scripts support generic centered 
spatial finite differencing, with optional upwinded differ- 
encing (and other off-center differencing) provided by ex- 
ternal macros. Currently we use second and fourth-order 
spatial finite differencing. 

The LazEv framework applies the MoL technique for 
solving first-order-in-time, hyperbolic PDE's, in which 
the PDE's are converted into coupled ODE's for all vari- 
ables at every gridpoint. This is achieved by choosing a 
discrete numerical grid and finite difference stencils for 
the spatial derivatives (the finite different stencils couple 
the fields at neighboring gridpoints), and then integrating 
the time derivatives of the fields at all gridpoints. This 
time integration can be carried out by standard ODE 
integrators. 

The LazEv MoL integrator provides a generic frame- 
work for integrating hyperbolic PDE's using Runge- 



Kutta or ICN style time integrations. The integra- 
tor itself has no knowledge of the system that is be- 
ing evolved. The MoL integrator provides internal 
timcbins pre-ministep, ministep, dissipation, and 
post-ministep, which are called during each step of 
the Runge-Kutta or ICN integration (there are several 
ministeps in each full timestep). The evolution system 
is chosen by registering routines with MoL to be evalu- 
ated during each of these bins. For example, in the BSSN 
system (see Sec. 1111 A0 we register routines to calculate 
the time derivative dt of all the BSSN variables in the 
ministep timebin, a routine to rcscale the BSSN confor- 
mal metric 7^ to unit determinant in post-ministep, 
and a routine to subtract off the numerical trace of the 
BSSN trace-free, conformal extrinsic curvature Ay in 
post-ministep. In the ADM system, on the other hand, 
we register routines only in the ministep bin. Additional 
evolution systems may also register routines during these 
steps. For example, the gauge evolution systems, which 
are independent of the main evolution systems, also reg- 
ister routines in ministep and post-ministep. 

Currently we have implemented ICN (second-order) 
and Runge-Kutta (second, third, and fourth-order). In 
some cases we use Kreiss-Oliger |45| dissipation of the 
form 

d t F -> RHS + (-l) r/2 e V +1 A+ r/2+1 A- r/2+X F, 

(2.1) 

where dtF = RHS is one of the evolution equations, hi 
is the gridspacing in the ith direction, Di + and are 
the forward and backwards differencing operators (in the 
ith direction) , and r is the order of the finite differencing 
used to evaluate RHS. 

We use several different styles of finite difference sten- 
cils, depending on the order of accuracy desired and the 
type of derivative considered. Our standard choices for 
fourth-order accurate evolutions are: 

+ 8F i+ ij : k — F i+2 j,k), (2.2) 

— 30F idt k + 16-Fi-i,,-, k - Fi_ 2<j , k), (2-3) 
9 ** Fi * h = lUdxdy 

[Fi-2,j-2,k — &Fi-l,j-2,k 
+ 8Fj + ij_2,fc — Fi + 2 : j-2,k 

— &{Fi-2,j—l,k — 8^i-l,j-l,fc 
+ 8-Fi+l,J-l,fc — ^i+2,j-l,fc) 
+ 8(-Fi-2,i+l,fe — ^Fi-l,j+l,k 

+ 8Fj + i j+i,fc — Fj-|_2,i+i,fe) 

— (Fi-2,j+2,k — 8F'i_i J+ 2 i fe 

+ 8F i+ ij +2 ,k — F i+2 ,j+2,k)] (2.4) 
Note that the mixed xy derivative is obtained by applying 
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the x-derivative and y-derivatives sequentially (the order 
is irrelevant). 

We do not use standard centered differencing for the 
advection terms (i.e. terms of form /3 l diF). For these 
terms we use the following upwinded stencils: 



d x Fi 



1 



i,j,k 



d x F, 



\2dx 

{—Fi-3,j,k + ~ 18-Pi-l,i,fc 

+ 10F i:j>k + 3F i+1 , jtk ) for »" • (2.5) 
1 



I2dx 

- 10F i:j>k - 3Fi- 1Jtk ) for [3 X > (2.6) 

In addition we use lower-order centered differencing 
on the planes adjacent to the boundaries. We have two 
choices for how the lower-order stencils are constructed. 
We can choose to use second-order centered differencing 
at all points on these planes and for all directions, or we 
can use second-order centered differencing only in the di- 
rection perpendicular to the boundary. When we use this 
latter choice we construct mixed second derivatives by 
applying the fourth-order accurate centered first deriva- 
tive operator (in a direction tangent to the boundary) to 
the second-order centered derivative operator (perpen- 
dicular to the boundary). 

Although we plan to implement excision boundary con- 
ditions soon, we currently use the puncture approach (see 
Ref. [42j for a comparison of methods), along with singu- 
larity avoiding slicings, to evolve black hole space-times. 
When using punctures, we found that we needed to either 
use second-order stencils (but still fourth-order Runge- 
Kutta time integration) in regions inside the apparent 
horizons, or use a second-order upwinded stencil for the 
advection terms over the entire computational domain. 
The former method produces significantly better wave- 
forms. See Sec. IV Cl for further details. 



III. FORMULATION 

Within our LazEv framework we have implemented 
support for several standard options in formulating nu- 
merical relativity modeling problems. Such a formula- 
tion includes selection of an Einstein evolution system 
and a choice of gauge for the evolving spacetime. In 
this section we discuss the formulation options presently 
available within our LazEv framework. Of these we have 
specific realizations which are appropriate for the fourth- 
order applications discussed in the subsequent sections. 



A. Evolution Systems 

Many alternative formulations of Einstein's evolution 
equations have been considered so far (see Ref. ^(| and 
references therein). Within the most popular Cauchy 



or 3 + 1 approach, one can currently choose among the 
first order symmetric hyperbolic formulations of the evo- 
lution equations which are mathematically very attrac- 
tive [l9l I20I I22I |23| . and various flavors of spatially 
second-order hyperbolic formulations which are numer- 
ically more tractable |47t l48j . While the former have the 
advantage of possessing a well-posed continuum limit in 
the presence of maximally dissipative boundary condi- 
tions, it is n ot y et clear how one can handle dozens of free 
parameters |2(| and dynamical gauge choices for the 
lapse and shift. On the other hand, some second-order 
hyperbolic formulations have proven to be empirically 
much more robust than others 0, 0, 0, • In this 
paper we will focus on the Baumgarte-Shapiro-Shibata- 
Nakamura (BSSN) n^.[T5l[T7| system, a strongly hyper- 
bolic system |5l|, |52| that has been shown to have some 
attractive stability properties [S3, E3 ■ 

The BSSN system is an extension of the standard 
ADM 54j system with better numerical properties than 
the original system. In the ADM system the Einstein 
Equations are split into six evolution equations for the 
metric (along with six auxiliary evolution equations for 
the extrinsic curvature) 



d Kij = 



= -2aK, 



13 > 

-DiDja, 



-a(Ri 



KK, 



2K lk K 



.))■ 



and four constraints equations 

H = R + K 2 - K i3 K ij = 0, 
M l = n,il\' J "''A.' = 0. 



(3.1) 
(3.2) 



(3.3) 
(3.4) 



Here do is the operator d t — £p, Cp is the Lie derivative 
with respect to the shift vector Di is the covariant 
derivative associated with the 3- metric 7^, Rij is the 
three-dimensional Ricci tensor, R the Ricci scalar, and 
K is the trace of Kij. 

The BSSN system of equations is obtained from the 
standard ADM equations by the following substitutions, 



7ij 
K i4 



^ A 



(3.5) 
(3.6) 



where 7 = det7ij = 1 and A\ = 0. Three additional 
variables 



r = -d 3 f ] 



(3.7) 



are also introduced. The BSSN variables (7^, A^-, K, <j>, 
and r l ) obey the following evolution equations (55| 



dojij = ~2aA 



13, 



do<j) = --^aK, 
6 

d A i:j = e-^ {-DiD ja + aR K 



(3.8) 
(3.9) 



■.TF 



4 



a (KAij - 2A ik A) y ) , (3.10) 
d K = -DWta + a (A^A^ + ^K 2 ^j , (3.11) 

d*p + Irdjft - 2A ij d ia + 
3 

2a (fijuA* + 6A^d j( p - \f 3 d 3 K^j , (3.12) 

where TF indicates that only the trace-free part of the 
tensor is used and Rij = Rij + R^ij is given by 

= -2D i D j ct>-2j ij D k D k <l> + 4:D i j ct>- 
4r/ ij D k <f>D k 4,, (3.13) 

Rij = -\l lm did m % + % (l d f) f k + f k f {ij)k + 

j lm (2T k l{i r j)km + f k im f , (3.14) 

and Di is the covariant derivative with respect to 7^ . T l 
is replaced by —djj^ in Eq's l|3.8|) - l|3.14|) wherever it 
is not differentiated. Note that Eq. I|3.12[l gives the dt 
derivative of f l rather than the do derivative. The Lie 
derivatives of the non-tensorial quantities (</>, 7^, and 
Ay) are given by 

= (3 k d k cj>+\d k f3 k , (3.15) 
6 

Cpjij = (3 k d k % + J ik djl3 k + 

l 3 kd^ k - ^ l3 d k [3 k , (3.16) 
CpAii = p k d k A l3 +A lk d 3 (3 k + 

A 3k d t f3 k - ^A l3 d k (3 k . (3.17) 

In addition to the evolution equations, the BSSN vari- 
ables must also obey the following seven differential con- 
straint equations 

U = R — AijA ij + \k 2 = 0, (3.18) 

o 

e^M 1 = djA ij + f* jk A? k + 6A ij dj<f> - 

~7 ij djK = 0, (3.19) 

g l = r + djf j = o, (3.20) 

and the following two algebraic constraint equations 

7 = 1, (3.21) 
A\ = 0. (3.22) 

We monitor the differential constraints but do not enforce 
them. However, we enforce the algebraic constraints by 
rescaling the evolved %j and subtracting off the trace of 
the evolved Aij at every timestep. 



We can evaluate (Rij) TF — Rij — l/3fijR in two ways. 
We could either calculate R^ and find its trace R nu- 
merically and subtract off the numerical trace, or we can 
use the Hamiltonian constraint Eq. (|3.18l) to calculate R 
and subtract that from Rij. Our tests showed very lit- 
tle difference between the two approaches even when the 
Hamiltonian constraint violations were relatively large. 
Unless otherwise specified we use the latter method. 



B. Gauge Choices 

We have implemented three basic types of lapse con- 
ditions: (i) maximal slicing, (ii) "Bona-Masso" type 

lapses [25[ |55l |5(| IE?} , (ii) and densitized lapses. The 
maximal slicing condition has the form 

Aa = (3 i d l K + aK l3 K' j , (3.23) 

and is implemented usin g a modified version of the 
Cactus 'Maximal' thorn |44| in combination with a 
fourth-order accurate version of Bernd Briigmann's 
'BAM_Elliptic' @ HI] thorn. This fourth-order accu- 
rate version of 'BAM_Elliptic' was developed at UTB by 
Mark Hannam. The "K-driver" lapses have the form 

d t a = ~a 2 f(a)A, (3.24) 

where 

d t A = d t K - £A, (3.25) 

and £ is some specified function (usually zero). We have 
also implemented modifications to this general form by 
replacing dta with doa, adding Di[3 l to the right-hand- 
side of Eq. (|3.24|) . and multiplying A by factors of -0" 
in Eq. (|3.24|) . where -0 is the conformal factor of the 
puncture data. Equation 13.2411 includes harmonic slic- 
ing (/(a) = 1), 1+log slicing (/(a) = 2/a), and 'shock- 
avoiding' slicing (/(a) = |(a(3 — a)) -1 ) |59l |. 

We have implemented 'Gamma-driver' shift condi- 
tions as well as static corotating shifts. The two ver- 
sions of the 'Gamma-driver' shifts that we implemented 
are 

= FB l -CP\ 
d t B l = d f r- V B\ (3.26) 

and 

d t p l = B l ~C{3\ 

d t B l = Fd t f l ~r]B l , (3.27) 

where 77 and £ are some prescribed functions over space, 
and F is a function of a and ip. Unless otherwise speci- 
fied, we take ( = 0. 
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C. Boundaries 

Thus far we have implemented fairly simple boundary 
conditions. For fourth-order upwinded stencils we alter- 
natively use centered finite differencing or second-order 
upwinding at the second point from the boundary, and 
second-order centered differencing at the first point from 
the boundary. The boundary points are filled using radia- 
tive boundary conditions [55( . However, when the shift is 
zero, or sufficiently small, we evolve 7^ on the boundary 
using the zero-shift form of Eq. (|3.8|) . We have observed 
that the boundary algorithm has very little effect on the 
convergence of the £ = 2 component of the waveforms in 
Sec. E| 

These boundary conditions are not known to be well 
posed and do lead to incoming constraint violating 
modes. Our future work will involve imposing constraint- 
preserving boundary conditions |48l l6Cl |6l| . 



IV. APPLES WITH APPLES TESTS 

We apply the LazEv framework to a standard set of 
numerical relativity tests known as 'Apples with Apples' 
tests |4l|. These tests consists of evolving spacetimes 
with B x T 3 topology with the advantage that there are 
no boundaries. The results from these testbeds indicate 
that LazEv is stable and fourth-order convergent with our 
choices of fourth-order stencils. 



A. Gauge wave test 

For this test we evolve the metric 

ds 2 = -H dt 2 + Hdx 2 + dy 2 + dz 2 , 

where 

,'2ir(x-tY 
H = 1 + A sin ' v ; 



(4.1) 



and d is the wavelength. The ADM variables for this 
metric have the form 



^xx — 

Ivy = 



1 + A sin 

Izz = 1, 

Att 



COS 



2ir{x - t) 
d 

/ 2ir(x-t) 
\ d 



1 + A sin 



2-n(x-t) 
d 



= i 1 + A sin 



2n(x - 1) 



(4.2) 
(4.3) 

(4.4) 
(4.5) 



where all remaining ADM variables (including (3" 1 ) are 
zero. Note that a obeys the harmonic slicing condition 
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FIG. 1: The L 2 norm of 8 lxx = -y xx - -y xx alytic , reseated 
by p 4 /16, for the one-dimensional 'Gauge Wave' test with 
A = 0.01. Note the near perfect overlap for 630 crossing 
times and that the larger p runs are convergent longer. The 
runs are acceptably accurate when the (non-rescaled) norm 
of the error is smaller than 10" 4 (i.e. 1% of A). 



In our simulations a is 'live' and we evolve it using the 
harmonic slicing condition. A two-dimensional test is 
obtained by the coordinate transformation 



y 



T2 {v ~ x) > 



Kxxj and Kyy 



d t a 



-a z K. 



(4.6) 



producing non-trivial ^ X x, 

In order to obtain stable runs we needed to include 
dissipation of the form given in Eq. (|2.1(l , with dissipation 
coefficient e = 0.125. 

These one and two-dimensional problems were evolved 
using a three-dimensional grid with periodic boundary 
conditions in all directions. We chose a wavelength (d) 
of 1 and constructed the grid so that it contained 6 + 1/h 
points (6 points for ghost-zones) in the non-trivial direc- 
tions and 9 points in the trivial directions (the stencil 
requires 3 ghost-zones), where h is the gridsize and 1/h 
is an integer. 

Our first test is a weak gauge wave with amplitude 
A = 0.01 and resolutions h = 0.02/p, where p = 2, 4, 8. 
Figure^shows the L2 norm of the error in ^ xx versus time 
(rescaled by a factor of p 4 /16). Note that fourth-order 
convergence (as evident by the overlap of the rescaled 
curves) implies that the error for the p = 8 case is 256 
times smaller than the p = 2 case. This relationship 
breaks down near 630 crossing times, though the higher 
resolution runs remain in the convergence regime longer. 
Using the criterion that the numerical solution is suffi- 
ciently accurate if the norm of the error is smaller than 
1% of the amplitude A, we find that the fourth-order 
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1000 



t (crossing times) 

FIG. 2: The L2 norm of Ti, rescaled by p 4 /16, for the one- 
dimensional 'Gauge Wave' test with A = 0.01. Note the near 
perfect overlap for 800 crossing times and that the larger p 
runs are convergent longer. The runs are acceptably accurate 
when the (non-rescaled) norm of the Hamiltonian constraint 
is smaller than 10~ 4 (i.e. 1% of A) 



code produces accurate results to t = 308, t = 720, and 
t = 865 for the p — 2, p — 4, and p — 8 runs respectively. 
Figure [21 shows the L2 norm of H (rescaled by a factor of 
p /16). Again, the higher resolution runs remain in the 
convergence regime longer. Using the criterion that the 
numerical solution is sufficiently accurate if the norm of 
the Hamiltonian constraint is smaller than 1% of the am- 
plitude A, we find that the fourth-order code produces 
accurate results to t = 130, t = 495, and t — 638 for 
the p — 2, p = 4, and p = 8 runs respectively. The 
Hamiltonian constrain violation is a stricter measure of 
the quality of the results than the error in j xx because 
the Hamiltonian constraint involves second derivatives of 
the metric (which are harder to model). 

Our next test is a stronger gauge wave. Figure shows 
the I/2 norm of Ti versus time (rescaled by a factor of 
p 4 /16) for the gauge- wave test with amplitude A = 0.1 
and resolutions h — 0.02/ p, where p = 2, 4, 8. Note 
that in this higher amplitude case the runs remain con- 
vergent for 80 crossing times (approximately one-tenth 
as long as the A = 0.01 runs). Using the criterion that 
the numerical solution is sufficiently accurate if the norm 
of the Hamiltonian constraint is smaller than 1% of the 
amplitude A, we find that the fourth-order code produces 
accurate results to t = 19.5, t = 46, and t — 60 for the 
p = 2, p — 4, and p = 8 runs respectively. Hence these 
A = 0.1 runs are accurate for roughly one-tenth the time 
that the A = 0.01 are accurate. 

Though we have applied our fully three-dimensional 
code to the last two test problems, the dynamics in these 
cases are non-trivial only in one (numerical grid) di- 
rection. Our last gauge wave test is non-trivial in two 
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FIG. 3: The L2 norm of TL, rescaled by p 4 /16, for the one- 
dimensional 'Gauge Wave' test with A = 0.1. Note the near 
perfect overlap for 80 crossing times and that the larger p 
runs are convergent longer. The runs are acceptably accurate 
when the (non-rescaled) norm of the Hamiltonian constraint 
is smaller than 10" 3 (i.e. 1% of A). 
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FIG. 4: The L2 norm of Ti., rescaled by p 4 /16, for the two- 
dimensional 'Gauge Wave' test with A = 0.1. Note the near 
perfect overlap for 55 crossing times and that the larger p 
runs are convergent longer. The runs are acceptably accurate 
when the (non-rescaled) norm of the Hamiltonian constraint 
is smaller than 10" 3 (i.e. 1% of A). 



grid directions. Figure 01 shows the L2 norm of Ti, ver- 
sus time (rescaled by a factor of p 4 /l6) for the two- 
dimensional gauge- wave test with amplitude A = 0.1 
and resolutions h = 0.02/ p, where p = 2, 4, 8. These 
two-dimensional runs are less stable than the correspond- 
ing one-dimensional run, crashing before 100 crossing 
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times. Convergence is, nonetheless, maintained for about 
55 crossing times, longer for the higher resolution runs. 
Using the criterion that the numerical solution is suf- 
ficiently accurate if the norm of the Hamiltonian con- 
straint is smaller than 1% of the amplitude A, we find 
that the fourth-order code produces accurate results to 
t = 10, t = 30, and t = 40 for the p = 2, p = 4, and 
p = 8 runs respectively. Note that these two-dimensional 
A = 0.1 runs are accurate roughly two-thirds as long as 
the corresponding one-dimensional A = 0.1 runs. 



B. Gowdy wave test 

In this section we present results for the 'Polarized 
Gowdy Wave' test. The 'Polarized Gowdy Wave' met- 
ric is given by 



ds 2 = t-V 2 e x / 2 (-dt 2 +dz 2 )+t(e p dx 2 +e- p dy 2 ) 
where 



P 
A 



Jo(27rf) C0S(27T2:), 

-27rfJ (27rf) Ji(27rf) cos 2 (2ttz) 
2tt 2 z 2 [j 2 (27rf) + Ji(27rf)] - 

±((2n) 2 {J 2 (2n) + J 2 (27r)}- 
27rJ (27r)Ji(27r)). 



(4.7) 



(4.8) 



(4.9) 



We use this metric to obtain initial data and evolve back- 
wards in time using the harmonic slicing condition. The 
above metric only varies as a function of z and t but 
we can obtain a two-dimensional problem by introducing 
new coordinates x — x, y — y, z = z + y. 

Here again we use a full three-dimensional grid to 
solve one and two-dimensional problems. We use peri- 
odic boundary conditions in all directions, and our grid 
consists of 1/h points (plus 6 for ghost-zones) in the 
non-trivial directions, and 9 points in the trivial direc- 
tions. Additionally, we needed to add dissipation of 
the form given in Eq. (|2.1|) . with dissipation coefficient 
e = 0.000625, to stabilize the runs. 

Figures |3] and show the fourth-order convergence of 
the Hamiltonian constraint and the ^ zz component of the 
metric of the one-dimensional Gowdy Wave test for 1000 
crossing times. Note that these convergence plots imply 
that the errors in j zz and Ti. in the p — 8 run are 256 
times smaller than these errors in the p — 2 run. Unlike 
the previous gauge- wave test, here the amplitude of the 
metric functions are damped. So our criterion for an 
acceptable value of the Hamiltonian constrain violation 
is time dependent. In this case we will use the criterion 
that the error norm, divided by the norm of j zz , must 
be smaller than 1%. At the end of the run the L2 of 
7 Z 2 (the function, not the error) is of order 1. Both the 
Hamiltonian constraint and the error in "f zz are smaller 
than this quality criterion for all resolutions. 

Figure [7| shows the fourth-order convergence of the 
Hamiltonian constraint for the two-dimensional Gowdy 
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FIG. 5: The L2 norm of Ti., rescaled by p 4 /16, for the one- 
dimensional 'Gowdy Wave' test. Note the good agreement 
between the curves for 1000 crossing times and that the evo- 
lution is backwards in time 
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FIG. 6: The L2 norm of the error in y zz , rescaled by the L2 
norm of 7 ZZ and by p 4 /16, for the one-dimensional 'Gowdy 
Wave' test. Note the good agreement between the curves 
for 1000 crossing times and that the relative error in y zz is 
less than 2 • 10" 4 at 1000 crossing times. The evolution is 
backwards in time. 



Wave test. The early time lack of convergence is due to 
roundoff effects, and, although the p = 2 curve does not 
lie on the p = 4 and p = 8 curves, the overlap of the p = 4 
and p — 8 curves confirm that the code is fourth-order 
convergent with sufficient resolution. 

Figure [S] shows the fourth-order convergence of ^f vv for 
the two-dimensional Gowdy Wave test. Unlike Ti, ^yy 
clearly demonstrates fourth-order convergence at the low 
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FIG. 7: The L2 norm of H, rescaled by p 4 /16, for the two- 
dimensional 'Gowdy Wave' test. The test was limited to 50 
crossing times due to limited resources. The early time lack of 
convergence is due to low amplitude high frequency noise in 
the p — 8 results. Although the two higher resolution curves 
do not lie on top of p — 2 curve, the overlap of the two high 
resolution curves indicates that the Hamiltonian constraint is 
converging to fourth-order. Note that the evolution is back- 
wards in time. 

(p = 2) resolution. 

We conclude from the gauge wave and Gowdy wave 
results that our BSSN code is stable, accurate, and con- 
vergent. 

V. HEAD-ON BINARY BLACK-HOLE 
COLLISIONS 

In this section we present results for head-on colli- 
sions of two equal-mass Misner-Wheeler-Brill-Lindquist 
(MWBL) black holes These results extend our 

previous tests of LazEv to more interesting nonlinear 
spacetimes containing binary black holes. 

The MWBL data represent conformally flat slices of 
multiple black hole space-times with n punctures. The 
data are parametrized by n puncture masses m, and n 
puncture positions ?i (in the conformal space). The ADM 
mass is given by the sum of m,, and the ADM linear 
momentum and angular momentum are zero. The data 
have the form 

7ij = if> 4 5ij, (5.1) 
Kij = 0, (5.2) 

where 

and Ti = \f — 



FIG. 8: The L2 norm of the error in j yy , rescaled by p 4 /16, for 
the two-dimensional 'Gowdy Wave' test. The test was limited 
to 50 crossing times due to limited resources. The excellent 
overlap of all three curves indicates that the code is converging 
to fourth-order. Note that the evolution is backwards in time. 

A. Setup 

We used MWBL octant-symmetric data consisting of 
two equal mass black holes aligned along the z-axis, al- 
lowing us to evolve the data using the {x > 0, y > 0, z > 
0} octant. In addition, we use a 'Transition Fisheye' 
transformation [lOj to enlarge the physical domain with- 
out sacrificing resolution near the punctures. The 'Tran- 
sition Fisheye' transformation is a smooth radial trans- 
formation from an inner resolution fixed by the Cactus 
gridspacing (h) to an outer resolution that is generally 
lower than the inner resolution. This transformation 
is parametrized by an outer 'de-resolution parameter' a 
(the effective grid-spacing at the edge of the grid is ah), 
a transition width parameter 's', and the center of the 
transition region £ ro'. The transformation has the form 

r physical = r(a + (1 - a)K{r)), (5.4) 

where r p ) lys i ca i is the physical radius corresponding to 
the coordinate radius r and lZ(r) is given by 

s /cosh^\ , 

*'>-*i^ ta UEW <5 ' 5) 

For these runs we set the mass parameter of the 
two holes to 0.5M, and the puncture positions to 
(0,0,-1.1515M) and (0, 0, 1.1515M). For most of the 
runs, the computational grid extended to 12.3M , which 
corresponds to a physical size of 26M. We also performed 
some short runs where the computational grid extended 
to 24. 6M, which corresponds to a physical size of 63M. 
These latter runs were used to determine the effect of 
the boundaries on the waveforms, constraint violations, 
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and horizon mass. The Fisheye parameters used in these 
runs were a — 3, s — 1.2M, r — 5.5M. We performed 
runs at resolutions of h = 1.1515M / (9p) with gridsizes of 
(96jo) 3 gridpoints along one octant, where p = 1,2, 4. The 
gridsizes were chosen so that the punctures lie halfway 
between gridpoints. For the runs with the boundary at 
63Af we used the same resolutions but chose gridsizes of 
(192p) 3 gridpoint along one octant, where p = 1,2. 

We evolved the data with the standard 1+log lapse, 
where the initial value of A (see Eq. 13.25(1 is zero, and a 
modified 1+log lapse, where the initial value of A is 

A(t = 0) = cexp(-(ip- iy 2 /a) . (5.6) 

We choose 1 for the initial value of the lapse. When using 
the modified lapse we set c = 0.5, and a — 0.8. This 
modified lapse, unlike the original 1+log lapse, collapses 
at the puncture in the continuum limit. Analytically the 
standard 1+log lapse should retain its initial value at 
the puncture throughout the entire evolution. Since we 
start with a lapse of one, the lapse should not collapse 
at the puncture. However, the lapse will collapse near 
the puncture, and a lack of resolution drives the lapse 
to zero at the puncture as well. However, as the grid is 
refined this artificial collapse at the puncture is delayed, 
leading to short wavelength features and blow-ups near 
the puncture. The initial value of A in Eq. (|5.6|l forces 
the lapse to collapse near the puncture in the continuum 
limit. Unfortunately this modification introduces its own 
problems as described below (see Fig. . 

We used the second form of the 'Gamma-Driver' shift 
Eq. (|3.27|l to evolve the shift. The initial value of the 
shift and its time derivative were zero. We used both a 
constant 77 (77 = 2.8/Af) and a spatially varying 77 equal 
to 2.8/M at infinity and 5.6/M at the punctures, as sug- 
gestedby Diener [64( • The spatially varying 77 was of the 
form 42] 



where r\ v and 77^ are parameters specifying 77 at the punc- 
ture and infinity respectively, and ip is the puncture data 
conformal factor. The function F in the 'Gamma-Driver' 
shift was set to F = |a/-0 4 . 

Pure fourth-order runs proved to be very unstable near 
the punctures both with and without upwinded stencils, 
and Kreiss-Oliger dissipation further destabilized these 
runs (near the punctures). We had success in stabilizing 
these runs with two techniques: (i) using second-order 
upwinding of advection terms over the entire grid, and (ii) 
a localized (within the apparent horizon) order reduction 
(LOR) to all second-order accurate spatial differentiation 
(including second-order upwinding of advection terms). 
This latter method provides significantly more accurate 
waveforms. When using the LOR method we found that 
lower-order accuracy is needed both at the punctures and 
at the origin. However, we only reduce the order of ac- 
curacy near the origin after a common apparent horizon 



forms (typically 5M after the common apparent horizon 
forms). We found that Kreiss-Oliger dissipation can be 
used to remove late time instabilities if the dissipation 
coefficient is set to zero in the LOR region. However, we 
did not use dissipation in the runs presented below. 

We implement LOR in the following way. First we 
calculate the time derivatives of all variables using the 
standard fourth-order stencils in Eq's 1(2. 2JI - 1(2.6(1 . Then 
we overwrite all points within some given coordinate- 
ellipsoid with the time derivatives calculated using the 
standard second-order stencils (with second-order up- 
winded advection terms). The dissipation operator in- 
side the LOR region can be either the same one chosen 
for the rest of the grid or a lower-order operator (in ei- 
ther case the dissipation coefficient inside the LOR region 
can be different from the dissipation coefficient used for 
the rest of the grid). Only the spatial finite differenc- 
ing is changed; the time integrator is the same for all 
gridpoints. The size of the ellipsoid is chosen at runtime 
and up to eleven ellipsoids may be used. The user is free 
to choose when each ellipsoid is activated. For example, 
in the head-on binary black hole collision case, we spec- 
ify a small ellipsoid for each of the individual apparent 
horizons and a larger ellipsoid for the common apparent 
horizon, where the larger ellipsoid is activated only after 
the common apparent horizons forms. The sizes of these 
ellipsoids can be determined by evolving with second- 
order accuracy over the entire grid and finding the sizes 
of the apparent horizons versus time. 

For the head-on-collision runs the LOR regions con- 
sisted of spheres of radius 0.512M centered on the punc- 
tures activated at t = 2M, as well as an ellipsoid 
with semi-axes {0.512M, 0.512M, 1.66M} activated at 
t = 11.6M (the punctures were located on the z-axis). 
For comparison, the common apparent horizon has a min- 
imum radius of 1.56M and a maximum radius of 1.88M 
at t = 11AM. The individual apparent horizons are ap- 
proximately spherical with radii 0.61A/ at t = 1.98M. 

We used radiative boundary conditions for all variables 
except (3 l (which were evolved via dtf3 l = B l on the 
boundary) . Table [I] gives the radiative boundary con- 
dition parameters for all variables. In Table |U 'a' is the 
Fisheye de-resolution parameter. We enforce the alge- 
braic constraints Eq's ((3.21(1 . 1(3.22(1 prior to applying the 
boundary conditions. Hence, the boundary points do not 
satisfy these constraints exactly. 

We updated the Zorro thorn of the Lazarus Toolkit 01 
to compute the Weyl scalars to fourth-order accuracy and 
to make it compatible with octant, bitant, and quadrant 
symmetries in any direction as well as tt rotation sym- 
metry about the z-axis. In addition we added a spherical 
harmonic decomposition routine to generate the 1 — 2 
and i = 4 modes shown below. 
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TABLE I: Radiative Boundary Condition Parameters 



Variables 



Asymptote 



speed 



9u 
Aij 

r 

K 
a 
B' 
A 



1 




(l/2)ln(a) 











1/a 

1/a 

1/a 

1/a 

y/2/a 

V2/a 

V2/a 

1/a 

1/a 



B. Second-Order Accurate Results 

We first confirmed that our code can reproduce the 
second-order accurate waveforms published in |Hfij . For 
these runs we used standard 1+log lapse (i.e. A(Q) = 0) 
and Gamma-Driver shift with 77 = 2.8. We ran with 
gridsizes of (96p) 3 gridpoints and grid resolutions of h = 
1.1515/(9p) for p = 1,2, 4. We calculated ip4 as well 
as its (I = 2 , to = 0) and {£ = 4, m = 0) components 
using Zorro [l0|. Note that the {£ = 2, to = 0) and 
{£ = 4, m = 0) modes of -04 are purely real for this test. 

Figure shows the (I = 2, to = 0) mode of ^4 at r = 
5M for these three resolutions along with the Richardson 
extrapolation of these data. In addition Fig. shows the 
differences 4('04| p= 2 — 04^=4) and ^4^=1 — ip4 1 p= 2- These 
two difference curves overlap reasonably well indicating 
that the waveforms are second-order accurate. However, 
the phase drift between the two curves makes an evalua- 
tion of the exact order of convergence difficult. Figure ITTjl 
shows the convergence rate for this mode given by 



= log; 



-0 4 |p=l ~ 4 , i\p=2 
04^=2 - "04|p=4 



The measured convergence rate oscillates wildly but av- 
erages to about 2. 

Strictly speaking r — 5-M is not in the radiation zone so 
■04 given above does not represent the asymptotic wave- 
form. However, the code's performance in calculating the 
[l = 2, to = 0) mode of ip4 at r = 5M is indicative of its 
performance in calculating this mode at large r. Addi- 
tionally, boundary errors contaminate the £ = 4 modes 
at large r (see Sec. IV Cj) . Thus by extracting at r — 5M 
we can maximize the amount of non-trivial, physically 
correct quasinormal oscillations in the waveforms. 

Figure ITU shows the Li norm of the Hamiltonian con- 
straint for these three runs. Note that the constraints 
have not been rescaled. From the figure we can see that 
the constraints tend to get bigger as the resolution is 
increased. The source of these constraint violations is lo- 
cated between the puncture and origin. High frequency 
features develop there, leading to strong constraint viola- 
tions. However, as seen in Fig. ^| these large constraint 
violations do not leak out of the apparent horizon. Thus, 
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FIG. 9: The top plot shows the {£ = 2, to = 0) mode of 
V>4 at r = 5M for the second-order evolution of the MWBL 
data for gridsizes h — 1.1515/(9p), as well as the Richard- 
son extrapolated value. The bottom plot shows the differ- 
ences (04|p=i —041^=2) and 4(^4 | p= 2 —^4|p=4) between these 
waveforms for this mode. Note that the latter difference has 
been rescaled by a factor of four. Second-order convergence 
is demonstrated by the reasonable overlap between these two 
differences. 



although there are significant problems inside the hori- 
zon, the region outside the horizon remains uncontami- 
nated. 

Our calculations of the gravitational waveforms from 
the Newman-Penrose scalar ip4, which contains higher- 
order spatial derivatives of the metric, do not reveal any 
additional degree of noise and distortion with respect 
the waveforms obtained from the Zcrilli-Moncricf formal- 
ism [34T.l42| . 
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FIG. 10: Convergence rate of the {I = 2, m = 0) mode of 
^4 at r= 5M for a purely second-order evolution of MWBL 
head-on collision data. 



0.6 



0.4 



0.2 











A = 

- p = 


i 

2 
4 








i \ 
i \ 

i \ 










/ 

.-■*' 








/ 










/ 
/ 









20 40 60 80 100 

t/M 



FIG. 11: The L2 norm of the Hamiltonian constraint viola- 
tion versus time for the second-order accurate head-on colli- 
sion runs. Note that the constraint violation increases with 
resolution. 



C. Fourth-Order Accurate Results 

Purely fourth-order runs of puncture data proved to be 
unstable with more than one black hole. We found that 
we could stabilize the evolution by reducing the spatial 
discretization to second-order inside the apparent hori- 
zons, and we successfully ran the MWBL head-on col- 
lision data with fourth-order discretization, fourth-order 
Runge Kutta time integration, and second-order LOR 
regions inside the apparent horizons. For these runs we 
used the spatially varying 77 form of the Gamma-driver 
shift given by Eq. (|3.27|) and Eq. (|5.7|l . FigurelP3ldemon- 
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FIG. 12: Hamiltonian constraint violation (absolute value) 
versus z along the z-axis at time t = 80M for the p = 4 second- 
order head-on collision run. Note that the large constraint 
violation inside the apparent horizon (r = 2.8M) has not 
leaked out into the exterior spacetime. 



strates the convergence of the (£ = 2, m = 0) component 
of V>4 at r = 5M. Note that the system is not strictly 
fourth-order accurate. There appears to be an additional 
phase discrepancy between the two differences. However, 
the amplitude of the differences appears to be falling at 
a rate consistent with better than fourth-order accuracy. 
The p — 4 run crashed at A7.5M due to an instability 
near the origin (see Fig. I19|l . This unstable mode was 
triggered by advection terms (the collapse of the lapse 
near the origin ensures that these are the only terms 
which can lead to a blowup). The fields which blew up 
most strongly were the T l . This blow-up does not appear 
to be related to the quadratic blow-up in <fr at the punc- 
ture discussed later. That term led to a blow up propor- 
tional to e (in TL) near the puncture, which masked the 
unstable behavior near the origin. Both blow-ups can be 
modified by various choices in the gauge conditions and 
it is likely that a better choice of gauge will lead to longer 
fourth-order evolutions. 

A fit of the (£ — 2, to = 0) mode of -04 to the 
quasinormal form / ~ e~ at sin(o>i) gave an exponen- 
tial damping factor of a = 0.084/M and frequency of 
uj — 0.373/Ajf. The exponential damping factor agrees to 
within 6% of the expected value of 0.0889625/Af, and the 
frequency agrees to within 0.2% of the expected value of 
0.373672/M 

In Fig. ^] we show the Richardson extrapolated value 
of the (i = 2, m = 0) mode of calculated using 
the second-order accurate results, along with the values 
obtained from second and fourth-order evolutions with 
LOR. Note that the medium resolution fourth-order run 
outperforms the high resolution second-order run. Fig- 
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FIG. 13: The top plot shows the (t = 2, m = 0) mode of ip4 at 
r = 5M, produced using fourth-order evolution with LOR, for 
3 resolutions, along with the Richardson extrapolation curve 
of these data. The p — 2 and p = 4 curves are indistinguish- 
able from the extrapolated curve. The bottom plot shows 
rescaled differences ip4,\ P =i — tp4\p=2 and 16(^4j P =2 — ipi\ p =A) 
(note the factor of 16). The convergence rate for the fourth- 
order runs using LOR inside the apparent horizons is not 
strictly fourth-order due to the phase discrepancy between the 
two curves, but average waveform differences between resolu- 
tions is consistent with fourth-order accuracy. Note that the 
amplitude of the dotted curve is less than the amplitude of 
the solid curve indicating better than fourth-order reduction 
in the amplitude of the errors. 



ure!15l shows a magnified view of Fig. ^] (note that the 
legends in the two figures are different). 

Figure shows the convergence of the (£ = 4, m = 
0) mode of ip4 at r = 5M calculated using fourth- 
order accuracy with LOR. The plot shows the differences 
"04 1 p= 1 — "04 1 p=2 and 16(V>4|p=2 — ip4\ P =i), as well the mode 
itself for p = 4. Note that the latter rescaled differ- 
ence is smaller than the former. This indicates that the 
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FIG. 14: A comparison of the {I = 2, m = 0) mode of ip4 
at r = 5M for second and fourth-order evolutions with LOR. 
Note that the p — 1 fourth-order waveform is better than the 
p = 2 second-order waveform. Similarly the p — 2 fourth- 
order waveform is better than the p = 4 second-order wave- 
forms. The Richardson extrapolated waveform (based on the 
second-order waveforms) is indistinguishable from the p = 2 
and p = 4 fourth-order waveforms. 



{£ = 4, m = 0) component of the waveform is converging 
faster than fourth-order. In addition, note that the error 
in the p = 1 waveform is comparable to the amplitude of 
the waveform. 

Figure IT71 shows a magnified view of the (£ = 4, m = 
0) mode of 1/14 at r = 5M for the p = 1, p = 2, and 
p = 4 fourth-order runs as well as the p = 4 second-order 
run. Note that the p — 4 second-order results are again 
inferior to the p = 2 fourth-order (with LOR) results. 

The £ = 4 modes, unlike the I — 2 modes, contain sig- 
nificant contamination from the boundary at late times. 
For example, when extracting at r — 13M we find that 
the amplitude of the (£ = 4, m = 4) mode of ip4 (which is 
non-zero due to boundary errors) is approximately 50% 
that of the (£ = 4,m = 0) mode at t = 45M. This con- 
tamination is stronger when the extraction sphere is fur- 
ther out. However, the waveform extracted at 5M shows 
significantly less contamination. The reason for this is 
that the incoming error is delayed by about 8M while the 
outgoing mode is advanced by 8M (for a net gain of 16M 
in reliable data). In addition the outgoing mode, which 
has a 1/r falloff, is correspondingly larger at this smaller 
radius. This boundary contamination problem can be 
mitigated by pushing the outer boundary further out. 
This can be achieved within the existing LazEv frame- 
work using a stronger fisheye transformation or adding 
more gridpoints. 

Although the waveforms from the fourth-order runs 
converge as expected the Hamiltonian constraint does 
not. High frequency features near the puncture and ori- 
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FIG. 15: A comparison of the (£ = 2, m = 0) mode of ^4 
at r = 5Af for second and fourth-order evolutions with LOR. 
Note that the p — 2 fourth-order waveform is better than the 
p = 4 second-order waveforms. The second-order p = 4 curve 
lies above the extrapolated curve while the fourth-order p — 1 
curve lies below. 



gin lead to large Hamiltonian constraint violation. At 
the points surrounding the puncture these high frequency 
features induce a quadratic blow-up in <\> for the p = 4 
run. This quadratic behavior, which is confined to the 

nearest neighbor points to the puncture, leads to a blow- 

.2 

up of e in the Hamiltonian constraint due to terms pro- 
portional to e 4 ^. This blow-up in 7i is entirely localized 
to the puncture and does not affect points outside. In 
the following figures we demonstrate that the large con- 
straint violations inside the apparent horizon do not leak 
out and that the constraint violations outside the appar- 
ent horizon converge to zero when boundary effects are 
removed. 

Figure IT8l shows the Hamiltonian constraint along the 
z-axis (the points surrounding the puncture have been 
removed) at t = 47.2M for the second and fourth-order 
runs for p — 4. The Hamiltonian constraint inside the 
horizon is as much as 10 4 times larger for the fourth- 
order run. However, outside the horizon the constraint 
violations are very similar. Note that we do not expect 
the constraints to converge to zero in this case because 
our boundary conditions are not constraint preserving 
and boundary constraint violations have contaminated 
the solution at this time (the boundaries were at 26M 
in physical coordinates). The p = 4 fourth-order run 
with LOR crashed due to an instability near the origin. 
Figure IT§1 shows the unstable mode in H. 

In Fig. [20| we show the L2 norms of the Hamilto- 
nian constraint TL, momentum constraint ./VP, and the 
BSSN constraints Q % (note that the x and y compo- 
nents of the both momentum and BSSN constraints are 
equal). We restricted these norms to the region outside 



FIG. 16: The (£ = 4, m = 0) mode of i> 4 for p = 4 as well as 
the differences between the p = 1 and p = 2 waveforms and 
the p — 2 and p = 4 waveforms (the latter rescaled by 16). 
The convergence rate of the (£ = 4, m = 0) mode of ^4 at r = 
5Af for the fourth-order runs with LOR is better than fourth- 
order. Note that the error in the p = 1 waveform (as evident 
by the difference between the p = 1 and p = 2 waveforms) is 
approximately 50% of the amplitude of the waveform at 20M 
and larger than the waveform beyond 30M. 



rphysicai = 3M (the horizon is at r physical ~ 2M) and 
inside r physical = 26M, where r physica i is the physical ra- 
dius. The outer boundaries were located at 63 M in phys- 
ical coordinates. These runs required double the number 
of gridpoints (along each direction) as the standard runs 
for the same resolutions. The reasonable overlap of the 
p = 1 and rescaled p = 2 curves for 43M of evolution 
indicates that the L2 norms are approximately fourth- 
order convergent (the convergence is lost near t = 40M 
due to constraint violating modes propagating in from 
the boundary). In addition, the amplitude of the con- 
straint violations is ~ 10~ 5 , which is nine order of mag- 
nitude smaller than the maximum Hamiltonian violation 
inside the horizon. We conclude from this figure that the 
second-order errors introduced by the LOR differencing, 
as well as the extremely large Hamiltonian constraint vi- 
olation inside the apparent horizon, do not propagate 
outside of the apparent horizon. 

The extreme Hamiltonian violation near the puncture 
can be removed using the modified 1+log lapse. How- 
ever, this modification tends to destabilize the runs at 
late times. Figure CHI shows the L2 norm of the Hamil- 
tonian constraint violation for the p = 4 fourth-order 
run with LOR. The modified lapse produces a constraint 
violation smaller than 1 for 25M but also causes the 
run to crash at 26. 5M. The standard 1+log lapse pro- 
duces a "stable" evolution to 47. 5M (see above) but has 
a constraint violation proportional to e at the punc- 
ture. If short term evolutions, like the ones required for 
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FIG. 17: The (£ = 4, m = 0) mode of */> 4 at r = 5M for 
fourth-order evolutions with LOR as well as second-order evo- 
lutions. Note that the p — 1 fourth-order waveform has an 
order 100% error and that the p = 2 fourth-order waveform 
has a substantially smaller error than the p = 4 second-order 
waveform. 



the Lazarus techniques, are required, then the modified 
1+log is preferred. 

Figure[52shows the horizon mass versus time, as calcu- 
lated by the AHFinderDirect thorn [6(j , for the p = 2 and 
p = 4 runs with LOR. The plot also shows the horizon 
mass of a p = 2 run with double the standard number of 
gridpoints in each direction and outer boundary at 63A/ 
in physical coordinates. The relatively sharp increase in 
the horizon mass at t = 32M is due to contamination 
from the boundary. When the boundary is moved out 
to 63Af (from 26M) this increase is delayed by roughly 
At = 37M. Also note that the slope of the increase is re- 
duced when the boundary is further out. The measured 
convergence rate of the horizon mass (based on runs with 
the boundary at 26M) exceeded 3.8 for 30M of evolution. 

The horizon-mass versus time plot shows that the late 
time boundary contamination contains significant erro- 
neous gravitational radiation. To determine when the 
calculated gravitational waveforms are good approxima- 
tions to the true waveforms, we need to confirm that (i) 
the waveforms converge and do not change significantly 
when the resolution is further increased, (ii) the wave- 
forms do not change significantly when the boundary ef- 
fects are removed (which we achieved by causally discon- 
necting the boundaries from the observer), and (iii) the 
Hamiltonian constraint violation (inside the extraction 
region) is much smaller than the waveform amplitude. 
We have shown that the waveforms converge, and that 
the Hamiltonian constraint violation converges in the ex- 
traction region (see Fig. I2(J|1 . Here we show the effects 
of the boundary on the Hamiltonian constraint and the 
waveforms. Figure l2"3l shows Hamiltonian constraint vio- 
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FIG. 18: The absolute value of the Hamiltonian constraint 
along the z-axis at t = 47.2A/ for the second and fourth-order 
runs. The points near the puncture have been removed (the 
violation at the puncture was 10 117 ) from the fourth-order 
data. Note that the constraint violation is 10 4 times bigger 
inside the horizon for the 4th order run, and that the extreme 
constraint violation does not leak out of the horizon (located 
at z = 2.8M). 
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FIG. 19: The unstable mode in Ti for the p = 4 fourth-order 
runs. 



lation along the z axis at time t = 40. 8M for the standard 
p — 2 run and a run with twice the standard number of 
gridpoints. The boundary in this latter case was at 63M 
(note that the z-axis plots the fisheye coordinate and that 
z fisheye = 12. 3 M corresponds to the physical coordinate 
Zphysicai = 26M and z fisheye = 24.6M corresponds to 
Zphysicai = 63Af). The figure shows the Hamiltonian con- 
straint outside the apparent horizon. Note that at this 
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FIG. 20: The L 2 norm, restricted to the region 3M < 
r v hyaicai < 26M, of all constraint violations. The bound- 
aries were located at 63M. All norms have been multiplied 
by 10 s and the p — 2 norms have been rescaled by an addi- 
tional factor of 16. In each panel the solid curves correspond 
to p = 1 and the dashed curves to p = 2. The reasonable 
overlap of the dashed and solid curves indicate that the con- 
straint violations are converging to zero to fourth-order. The 
constraints no longer converge to zero after t ~ 43M due to 
boundary effects. 
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FIG. 22: The horizon mass versus time for the fourth-order 
runs with LOR. The thick solid line and the dotted line show 
the horizon mass versus time for the p = 2 resolution with the 
physical boundary at 63M and 26M respectively. Note that 
the sharp increase in mass at later times is due to boundary 
effects. The dashed line shows the horizon mass for the p — 4 
run with the boundaries at 26M. The error in the horizon 
mass at t — 63. 5M is 2.6% for the p — 2 run with boundary 
at 26M. 
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FIG. 21: The L2 norm of the Hamiltonian constraint viola- 
tion for the p = 4 fourth-order evolution with LOR for the 

.2 

standard and modified 1+log lapses. Note the e blow-up 
in the run using the standard 1+log lapse and that the run 
using the modified 1+log lapse crashes at 26. 5M. 



time the boundary errors for the larger run have just 
reached Zfi S heye = 12M. The Hamiltonian constraint, in 
the range Zf lsheye = 5M to z f . lsheye = 12M, is 100 times 
smaller for the run with the larger boundary. Even at 
t = 70-A/, when the boundary errors have contaminated 
the entire grid, the Hamiltonian constraint is 20 times 
smaller for the larger run in this range. 

As seen in Fig. |2J the effect of reducing the boundary 
noise is not readily apparent in the (£ = 2,m = 0) mode 
of ip4. However, as seen in Fig. the effect is readily 
observable in the (£ = 4, m = 0) mode. We conclude, 
based on these two figures, that the (£ = 2, m = 0) mode 
is represented accurately to 55M when the boundaries 
are located at 26M, while the (I — 4, m — 0) mode is 
represented accurately only to 33M for the same bound- 
ary location. The (£ = 2, m = 0) mode is more accurate 
because the boundary contamination contains relatively 
high frequencies which are filtered more effectively by the 
(I = 2, m = 0) angular integration. Note that we placed 
the observer at r — 5M and that the waveforms become 
inaccurate sooner when the extraction sphere is placed 
at larger radii. 

We conclude this section by showing the waveforms 
produced with fourth-order centered spatial finite dif- 
ferencing for all terms except the advection terms, for 
which we used second-order upwindcd differencing. The 
waveforms produced by this method are inferior to those 
produced with LOR techniques but superior to those pro- 
duced by the standard second-order technique. The time 
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FIG. 23: The Hamiltonian constraint along the z axis at t = 
40. 8M for the fourth-order (with LOR) runs with gridspac- 
ing h = 1.1515/18 (i.e. p = 2). Zf is h e ye = 12. 3M corresponds 
to 26M in physical coordinates and Zfi S he ye = 24. 6M corre- 
sponds to 63A/ in physical coordinates. Boundary contam- 
ination for the larger run has just reached Zfi S h ey e = VIM. 
Note that in the range 5M < zfi S h ey e < VIM the Hamil- 
tonian constraint is 100 time smaller when the boundary is 
pushed out to 63M. 
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FIG. 25: The (i = 4,m = 0) mode of V>4 (observer at r = 5M) 
for the fourth-order (with LOR) p = 2 run with the physical 
boundary at 26M (standard) and 63M. The effect of the 
boundary is significant at t > 33M. 
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the second-order accuracy of the algorithm, the wave- 
forms produced with this technique are more accurate 
than those produced by the purely second-order spatial 
differencing. Figure I77I shows the {£ = 2, m = Q) mode 
of -04 at r = 5M for the standard second-order evolution 
as well as the mixed fourth-order with second-order up- 
winding evolution. Note that the p — 2 waveform from 
the latter technique is of a similar quality to the p = 4 
waveform from the purely second-order technique. The 
p — 4 run crashed at 47. 5M with the same mode that 
killed the p = 4 LOR run. Both techniques crashed at 
the same time and in the same way because the insta- 
bility near the origin is driven by advection terms, and 
both techniques use the same upwinded advection sten- 
cil near the origin. Figure [nil shows the unstable mode 
(along the z-axis) that crashes the p = 4 fourth-order 
runs (with LOR and with second-order upwinding). 



FIG. 24: The (£ = 2,m = 0) mode of ip 4 (observer at r = 5M) 
for the fourth-order (with LOR) p — 2 run with the physical 
boundary at 26M (standard) and 63M. The effect of the 
boundary is not significant until t — 55M. 



integration was carried out with the standard fourth- 
order Runge Kutta integrator. Figure|2fj]shows the differ- 
ences ip4\ p= i — ip4:\p=2 and ip4,\p=2 — ip4\p=4 with the latter 
rescaled by 4 (to indicate second-order convergence) for 
the (£ = 2, m = 0) mode of ^4 at r = 5M. The two 
curves do not overlap exactly due to phase drift. Despite 



We computed the energy radiated from the head-on 
collision by several different methods. We first estimated 
it by computing the difference between the final hori- 
zon mass and the total ADM mass, and obtained a ra- 
diated energy in the range (7 — 8) x 10~ 4 Af. We also 
used the Lazarus method to extract Cauchy data for 
the Teukolsky equation at relatively early evolution times 
T < 20M, and obtained E radlated = (8 ± 1) x 1CT 4 A/. 
The direct integration of the I — 2 and £ = 4 wave- 
forms that we present here, extracted at numerical radii 
between 5M and 10M, produce energy estimates of the 
order ~ 6.6 x 10 _4 M. 
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FIG. 26: Convergence plot for the mixed fourth-ordered cen- 
tered differencing with second-order upwinded differencing 
runs. Note that the difference ip4,\ p =2 — ip4\ P =4 has been multi- 
plied by 4, indicating approximate second-order convergence. 
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FIG. 27: A comparison of the {I — 2, m = 0) mode of i/>4 
at r — 5Af produced using the standard second-order evolu- 
tion and a mixed fourth-order with second-order upwinding 
evolution. Note that the p — 2 waveform from the mixed 
fourth-order/second-order upwinding is of similar quality to 
the p = 4 second-order waveform. The extrapolated values 
used in this plot are based on the waveforms from the fourth- 
order runs with LOR. 



VI. CONCLUSION 

We developed a new framework, LazEv, for evolving 
the Einstein equations using 3+1 decompositions. LazEv 
is capable of evolving with arbitrary-order finite differ- 
ence stencils along with second, third, and fourth-order 
time integrators. The overall LazEv design has a few 



novel features which will improve significantly upon pre- 
vious setups which have been used in the Lazarus ap- 
proach 0, [13, EH E3, 

LazEv is a flexible and modu- 
lar evolution package, well-suited for rapid development 
and experimentation with well-posed hyperbolic systems 
of evolution equations, new 'live' gauge conditions, and 
sophisticated boundary conditions using arbitrary order 
finite differencing. 

We implemented the BSSN formulation using this new 
framework and demonstrated that this new code passes 
the Apples with Apples testsuite and that the code re- 
produces the second-order accurate head-on binary black 
hole collisions waveforms published in Ref. |55|, We 
found that fourth-order accurate evolutions of the same 
data were not stable without some reduction of order 
and that the most accurate waveforms were obtained 
by evolving the data with second-order accurate sten- 
cils inside the apparent horizons and fourth-order sten- 
cils outside the horizon. In that case we found that 
the waveforms were fourth-order convergent and that 
the fourth-order accurate 192 3 gridpoint runs outper- 
formed the second-order- accurate 384 3 gridpoint runs. 
This means that our fourth-order evolutions give bet- 
ter quality waveforms with over an order-of-magnitude 
smaller computational expense when compared to the 
second-order evolutions. 

We also found that using a second-order upwinded 
stencil for the advection terms was sufficient to stabilize 
the runs. However, this introduces a second-order error 
in the waveform. Nevertheless, the waveforms produced 
by this latter method are superior to those produced by 
ordinary second-order accurate evolutions. 

In this paper we demonstrated that the LazEv frame- 
work can be used to evolve binary black hole spacetimes. 
We plan to extend this work to include orbiting black 
holes starting from initial data based on the conformal 
thin-sandwich formulation [6^ | as well as study recoil ve- 
locities from unequal mass head-on collisions We 
also plan to extend the framework to include excision, 
fixed mesh refinement, constraint damping [7(1 l7lj|. and 
constraint preserving boundary conditions. 
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